#include "stdafx.h"
#include <cuda_runtime.h> 
#include <cutil_inline.h> 

//*****************************************************************************
__constant__ double		d_interHtc[NUM_MATERIAL * NUM_MATERIAL];
__constant__ double		d_density[NUM_MATERIAL];
__constant__ double		d_sh[NUM_MATERIAL];
dieNum2D				dieNum[NUM_DIE];
plateNum2D				plateNum;
int						airThickNum;
int						row, col;

//*****************************************************************************
int initInterHtc(material substance[], double dx, double *interHtc)
{
	double	r, k;
	for(int i = 0; i < NUM_MATERIAL; i++)
		for(int j = i; j < NUM_MATERIAL; j++)
		{
			if (substance[i].both[j] == 0) r = 0; else r = 1 / substance[i].both[j];
			if (i == 0 && (j == 1 || j == 2)) k = TOUCH_FACTOR; else k = 1;
			interHtc[i * NUM_MATERIAL + j] = interHtc[j * NUM_MATERIAL + i] = 
				k * dx / (r + dx / 2 / substance[i].htc + dx / 2 / substance[j].htc);
			if (i == 0 && (j == 1 || j == 2)) 
				interHtc[i * NUM_MATERIAL + j] = interHtc[j * NUM_MATERIAL + i] = 7;
		}
		return 0;
}

dieNum2D getDieNum(die2D die, double dx)
{
	dieNum2D dieNum;
	dieNum.heightNum = (int)(die.height /dx);
	dieNum.widthNum = (int)(die.width / dx);
	return dieNum;
}

plateNum2D getPlateNum(plate2D plate, double dx)
{
	plateNum2D plateNum;
	plateNum.heightNum = (int)(plate.height /dx);
	plateNum.widthNum = (int)(plate.width / dx);
	return plateNum;
}

//*****************************************************************************
int setConstant(material substance[], double dx)
{
	double		h_interHtc[NUM_MATERIAL * NUM_MATERIAL];	
	initInterHtc(substance, dx, h_interHtc);
	cudaMemcpyToSymbol(d_interHtc, h_interHtc, 
		sizeof(double) * NUM_MATERIAL * NUM_MATERIAL, 0, cudaMemcpyHostToDevice);

	double		h_density[NUM_MATERIAL], h_sh[NUM_MATERIAL];
	for(int i = 0; i < NUM_MATERIAL; i++)
	{
		h_density[i] = substance[i].density;
		h_sh[i] = substance[i].sh;
	}
	cudaMemcpyToSymbol(d_density, h_density, sizeof(double) * NUM_MATERIAL, 0, cudaMemcpyHostToDevice);
	cudaMemcpyToSymbol(d_sh, h_sh, sizeof(double) * NUM_MATERIAL, 0, cudaMemcpyHostToDevice);
	return 0;
}

//*****************************************************************************
int getGrid(int size)
{
	return (size + BLOCK_SIZE - 1) / BLOCK_SIZE;
}

//*****************************************************************************
__global__ void kernelInitBlock(int col, int x, int y, int width, int size,
								int substanceId, double initTemp, 
								int	*d_grid, double *d_tempCurr, double *d_tempNext)
{
	int id = blockIdx.x * blockDim.x + threadIdx.x;
	if (id < size)
	{
		int r = id / width;
		int c = id - r * width;
		int location = (x + r) * col + (y + c);
		d_grid[location] = substanceId;
		d_tempCurr[location] = initTemp;
		d_tempNext[location] = initTemp;
	}
}

int initGridCuda(material substance[], int *d_grid, double *d_tempCurr, double *d_tempNext)
{
	int x, y, length, width, k;
	dim3 block(BLOCK_SIZE);
	//init block - Gas - Air
	k = 3;
	x = 0; length = row; y = 0; width = col; 
	dim3 totalGrid(getGrid(length * width));
	kernelInitBlock<<<totalGrid, block>>>(col, x, y, width, length * width, 
		substance[k].id, substance[k].initTemp, d_grid, d_tempCurr, d_tempNext);

	//init block - Plate - 
	k = 0;
	x = airThickNum + dieNum[0].heightNum; length = plateNum.heightNum; 
	y = (col - plateNum.widthNum) / 2; width = plateNum.widthNum; 
	dim3 plateGrid(getGrid(length * width));
	kernelInitBlock<<<plateGrid, block>>>(col, x, y, width, length * width, 
		substance[k].id, substance[k].initTemp, d_grid, d_tempCurr, d_tempNext);

	//init block - Die Up -  
	k = 1;
	x = airThickNum; length = dieNum[0].heightNum; 
	y = (col - dieNum[0].widthNum) / 2; width = dieNum[0].widthNum;
	dim3 dieUpGrid(getGrid(length * width));
	kernelInitBlock<<<dieUpGrid, block>>>(col, x, y, width, length * width, 
		substance[k].id, substance[k].initTemp, d_grid, d_tempCurr, d_tempNext);

	//init block - Die Down - 
	k = 2;
	x = airThickNum + dieNum[0].heightNum + plateNum.heightNum; length = dieNum[1].heightNum; 
	y = (col -  dieNum[1].widthNum) / 2; width = dieNum[1].widthNum; 
	dim3 dieDownGrid(getGrid(length * width));
	kernelInitBlock<<<dieDownGrid, block>>>(col, x, y, width, length * width, 
		substance[k].id, substance[k].initTemp, d_grid, d_tempCurr, d_tempNext);


	return 0;
}

//*****************************************************************************
__global__ void kernelTimeStep(int col, int col_2, int matrixSize, double dxdx, 
							   int *d_grid, double *d_timeStep)
{
	int id = blockIdx.x * blockDim.x + threadIdx.x;
	if (id < matrixSize)
	{
		int r = id / col_2;
		int c = id - r * col_2;
		int location = (1 + r) * col + (1 + c);
		int k = d_grid[location];
		d_timeStep[id] = d_density[k] * d_sh[k] * dxdx
						/ (d_interHtc[k * NUM_MATERIAL + d_grid[location - 1]] 
						+ d_interHtc[k * NUM_MATERIAL + d_grid[location + 1]]
						+ d_interHtc[k * NUM_MATERIAL + d_grid[location - col]] 
						+ d_interHtc[k * NUM_MATERIAL + d_grid[location + col]]);
	}
}

//kernel of min of an Array
__global__ void kernelGetMin(const int count, double *d_array)
{
	int id = blockIdx.x * blockDim.x + threadIdx.x;
	if (id < count)
		if (d_array[id] > d_array[id + count]) d_array[id] = d_array[id + count];
}

double getArrayMin(double *d_array, int arrayLength)
{
	dim3		block(BLOCK_SIZE);
	dim3		grid(getGrid(arrayLength / 2));
	int			count = (int)(arrayLength / 2);
	while(count > 1)
	{
		kernelGetMin<<<grid, block>>>(count, d_array);
		count = (count + 1) / 2;
	}
	double		first, last;
	cudaMemcpy(&first, d_array, sizeof(double), cudaMemcpyDeviceToHost);
	cudaMemcpy(&last, d_array + arrayLength - 1, sizeof(double), cudaMemcpyDeviceToHost);
	
	return first < last ? first : last;
}

double getMaxTimeStep(double dx, int *d_grid)
{
	
	double		*d_timeStep;
	int			matrixSize = (row - 2) * (col - 2);
	cudaMalloc((void**)&d_timeStep, matrixSize * sizeof(double));
	dim3		block(BLOCK_SIZE);
	dim3		grid(getGrid(matrixSize));

	kernelTimeStep<<<grid, block>>>(col, col - 2, matrixSize, dx * dx, d_grid, d_timeStep);

	////****************************************
	//double		*h_test = (double *)malloc(sizeof(double) * matrixSize);
	//cudaMemcpy(h_test, d_timeStep, sizeof(double) * matrixSize, cudaMemcpyDeviceToHost);
	//output(h_test, row - 2, col - 2, 5, 0, 0);
	//free(h_test);
	////****************************************

	return getArrayMin(d_timeStep, matrixSize);
}

//*****************************************************************************
__global__ void kernelCalcNextTemp(double timeStep, int col, int col_2, int matrixSize, double dxdx,
								   int *d_grid, double *d_tempCurr, double *d_tempNext)
{
	int id = blockIdx.x * blockDim.x + threadIdx.x;	
	if (id < matrixSize)
	{
		int r = id / col_2;
		int c = id - r * col_2;
		int location = (1 + r) * col + (1 + c);
		int k = d_grid[location];
		double temp = d_tempCurr[location];
		d_tempNext[location] = d_tempCurr[location] + timeStep / d_density[k] / d_sh[k] / dxdx * 
			(d_interHtc[k * NUM_MATERIAL + d_grid[location - 1]] * (d_tempCurr[location - 1] - temp) +
			d_interHtc[k * NUM_MATERIAL + d_grid[location + 1]] * (d_tempCurr[location + 1] - temp) +
			d_interHtc[k * NUM_MATERIAL + d_grid[location - col]] * (d_tempCurr[location - col] - temp) +
			d_interHtc[k * NUM_MATERIAL + d_grid[location + col]] * (d_tempCurr[location + col] - temp));
	}
}

double calcTempFieldCuda(double timeStep, double dx, int *d_grid, double *d_tempCurr, double *d_tempNext)
{
	int			pointX = airThickNum + dieNum[0].heightNum + plateNum.heightNum / 2;
	int			pointY = col / 2;
	double		totalTime = 0;
	double		*tag;
	int			matrixSize = (row - 2) * (col - 2);
	dim3		block(BLOCK_SIZE);
	dim3		totalGrid(getGrid(matrixSize));

	//int			matrixTotalSize = row * col;
	//int			matrixMemSize = sizeof(double) * matrixTotalSize;
	//double		*h_tempCurr = (double *)malloc(matrixMemSize);
	//cudaMemcpy(h_tempCurr, d_tempCurr, matrixMemSize, cudaMemcpyDeviceToHost);
	//h_tempCurr[0] = PLATE_INIT_TEMP;
	//output(h_tempCurr, row, col, totalTime, pointX, pointY);

	double		tempGoal = GOAL_TEMPERATURE;
	double		keyPointTemp;
	cudaMemcpy(&keyPointTemp, d_tempCurr + pointX * col + pointY, sizeof(double), cudaMemcpyDeviceToHost);

	while(keyPointTemp > tempGoal)
	{
		totalTime += timeStep;

		kernelCalcNextTemp<<<totalGrid, block>>>(timeStep, col, col - 2, matrixSize, dx * dx, d_grid, d_tempCurr, d_tempNext);
		tag = d_tempCurr; d_tempCurr = d_tempNext; d_tempNext = tag;

		//if ((int)(totalTime / OUTPUT_FREQUENCY) > (int)((totalTime - timeStep) / OUTPUT_FREQUENCY)) 
		//{
		//	cudaMemcpy(h_tempCurr, d_tempCurr, matrixMemSize, cudaMemcpyDeviceToHost);
		//	h_tempCurr[0] = PLATE_INIT_TEMP;
		//	output(h_tempCurr, row, col, totalTime, pointX, pointY);
		//	keyPointTemp = h_tempCurr[pointX * col + pointY];
		//}
		//else
			cudaMemcpy(&keyPointTemp, d_tempCurr + pointX * col + pointY, sizeof(double), cudaMemcpyDeviceToHost);
		//printf("%15.3lf  %15.3lf\n", totalTime, keyPointTemp);
	}
	//cudaMemcpy(h_tempCurr, d_tempCurr, matrixMemSize, cudaMemcpyDeviceToHost);
	//h_tempCurr[0] = PLATE_INIT_TEMP;
	//output(h_tempCurr, row, col, totalTime, pointX, pointY);

	return totalTime;
}

//********************************************************************************
double getCostTime(LARGE_INTEGER start, LARGE_INTEGER end)
{
	LARGE_INTEGER nFrequency;
	QueryPerformanceFrequency(&nFrequency);
	double time = (double)TIME_ACCURACY * (end.QuadPart - start.QuadPart) / nFrequency.QuadPart;
	return time / TIME_ACCURACY;
}


//*****************************************************************************
__global__ void getConstant(double *test)
{
	int i;
	for(i = 0; i < NUM_MATERIAL * NUM_MATERIAL; i++)
		test[i] = d_interHtc[i];

	return;
}
__global__ void getConstant1(double *test)
{
	int i;
	for(i = 0; i < NUM_MATERIAL; i++)
		test[i] = d_density[i];

	return;
}
__global__ void getConstant2(double *test)
{
	int i;
	for(i = 0; i < NUM_MATERIAL; i++)
		test[i] = d_sh[i];

	return;
}
//*****************************************************************************
double fourierCuda2D(int argc, char** argv, 
						struct material substance[], struct die2D die[], 
						struct plate2D plate, double airThick, double dx)
{
	 //init cuda device
	if (cutCheckCmdLineFlag(argc, (const char**)argv, "device"))
		cutilDeviceInit(argc, argv);
	else
		cudaSetDevice(cutGetMaxGflopsDeviceId()); 
	//timer
	LARGE_INTEGER	nStartCounter, nStopCounter;
	double			nTime;
	//timer starts
	QueryPerformanceCounter(&nStartCounter);

	setConstant(substance, dx);

	////****************************************
	////test d_interHtc
	//double *d_testHtc;
	//cudaMalloc((void**) &d_testHtc, sizeof(double) * NUM_MATERIAL * NUM_MATERIAL);
	//getConstant<<<1, 1, 0>>>(d_testHtc);
	//double *h_testHtc = (double *)malloc(sizeof(double) * NUM_MATERIAL * NUM_MATERIAL);
	//cudaMemcpy(h_testHtc, d_testHtc, sizeof(double) * NUM_MATERIAL * NUM_MATERIAL, cudaMemcpyDeviceToHost);
	//for(int i = 0; i < NUM_MATERIAL * NUM_MATERIAL; i++)
	//	printf("interHtc @ device = %lf \n", h_testHtc[i]);
	//free(h_testHtc);
	//cudaFree(d_testHtc);
	////****************************************
	////test d_density
	//double *d_testDensity;
	//cudaMalloc((void**) &d_testDensity, sizeof(double) * NUM_MATERIAL);
	//getConstant1<<<1, 1, 0>>>(d_testDensity);
	//double *h_testDensity = (double *)malloc(sizeof(double) * NUM_MATERIAL);
	//cudaMemcpy(h_testDensity, d_testDensity, sizeof(double) * NUM_MATERIAL, cudaMemcpyDeviceToHost);
	//for(int i = 0; i < NUM_MATERIAL; i++)
	//	printf("density @ device = %lf \n", h_testDensity[i]);
	//free(h_testDensity);
	//cudaFree(d_testDensity);
	////****************************************
	////test d_sh
	//double *d_testSh;
	//cudaMalloc((void**) &d_testSh, sizeof(double) * NUM_MATERIAL);
	//getConstant2<<<1, 1, 0>>>(d_testSh);
	//double *h_testSh = (double *)malloc(sizeof(double) * NUM_MATERIAL);
	//cudaMemcpy(h_testSh, d_testSh, sizeof(double) * NUM_MATERIAL, cudaMemcpyDeviceToHost);
	//for(int i = 0; i < NUM_MATERIAL; i++)
	//	printf("sh @ device = %lf \n", h_testSh[i]);
	//free(h_testSh);
	//cudaFree(d_testSh);
	////****************************************

	dieNum[0] = getDieNum(die[0], dx);
	dieNum[1] = getDieNum(die[1], dx);
	plateNum = getPlateNum(plate, dx);
	airThickNum = (int)(airThick / dx);

	row = dieNum[0].heightNum + dieNum[1].heightNum + plateNum.heightNum + airThickNum * 2;
	col = (dieNum[0].widthNum > dieNum[1].widthNum ? dieNum[0].widthNum : dieNum[1].widthNum) + 2 * airThickNum;
	//printf("row = %d \t col = %d \n", row, col);

	int			*d_grid;
	double		*d_tempCurr;
	double		*d_tempNext;
	int			matrixSize = row * col;
	cudaMalloc((void**)&d_grid, matrixSize * sizeof(int));
	cudaMalloc((void**)&d_tempCurr, matrixSize * sizeof(double));
	cudaMalloc((void**)&d_tempNext, matrixSize * sizeof(double));

	initGridCuda(substance, d_grid, d_tempCurr, d_tempNext);
	////****************************************
	//int			*h_testInt = (int *)malloc(sizeof(int) * row * col);
	//double		*h_test = (double *)malloc(sizeof(double) * row * col);
	//cudaMemcpy(h_testInt, d_grid, sizeof(int) * row * col, cudaMemcpyDeviceToHost);
	//outputInt(h_testInt, row, col, 10, 0, 0);
	//cudaMemcpy(h_test, d_tempCurr, sizeof(double) * row * col, cudaMemcpyDeviceToHost);
	//output(h_test, row, col, 1, 11, 0);
	//cudaMemcpy(h_test, d_tempNext, sizeof(double) * row * col, cudaMemcpyDeviceToHost);
	//output(h_test, row, col, 2, 12, 0);
	//free(h_test);
	//free(h_testInt);
	////****************************************

	//time step
	double maxTimeStep = getMaxTimeStep(dx, d_grid);
	//printf("Max Time Step = %f \n", maxTimeStep);
	double timeStep = TIME_STEP_FACTOR * maxTimeStep;

	double totalTime = calcTempFieldCuda(timeStep, dx, d_grid, d_tempCurr, d_tempNext);

	//timer2 ends
	QueryPerformanceCounter(&nStopCounter);
	//get time
	nTime = getCostTime(nStartCounter, nStopCounter);

	printf("totalTime = %lf\n", totalTime);
	printf("totalStep = %d\n", (int)(totalTime / timeStep));
	printf("Compute Cycle = %f\n", nTime);

	return 0;
}